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ABSTRACT 

We study the properties of ISM substructure and turbulence in hydrodynamic (AMR) 
galaxy simulations with resolutions up to 0.8 pc and 5 x 10 3 Mq. We analyse the 
power spectrum of the density distribution, and various components of the velocity 
field. We show that the disk thickness is about the average Jeans scale length, and is 
mainly regulated by gravitational instabilities. From this scale of energy injection, a 
turbulence cascade towards small-scale is observed, with almost isotropic small-scale 
motions. On scales larger than the disk thickness, density waves are observed, but there 
is also a full range of substructures with chaotic and strongly non- isotropic gas velocity 
dispersions. The power spectrum of vorticity in an LMC-sized model suggests that an 
inverse cascade of turbulence might be present, although energy input over a wide 
range of scales in the coupled gaseous+stellar fluid could also explain this quasi-2D 
regime on scales larger than the disk scale height. Similar regimes of gas turbulence are 
also found in massive high-redshift disks with high gas fractions. Disk properties and 
ISM turbulence appear to be mainly regulated by gravitational processes, both on large 
scales and inside dense clouds. Star formation feedback is however essential to maintain 
the ISM in a steady state by balancing a systematic gas dissipation into dense and 
small clumps. Our galaxy simulations employ a thermal model based on a barotropic 
Equation of State (EoS) aimed at modelling the equilibrium of gas between various 
heating and cooling processes. Denser gas is typically colder in this approach, which 
is shown to correctly reproduce the density structures of a star-forming, turbulent, 
unstable and cloudy ISM down to scales of a few parsecs. 

Key words: galaxies: evolution - galaxies: ISM - galaxies: star formation - galaxies: 
structure - ISM: kinematics and dynamics - ISM: structure 



1 INTRODUCTION 

Observations of the interstellar medium (ISM) of nearby 
galaxies reveal a large variety of complex substructures, 
such as density waves, molecular clouds, filamentary struc- 
tures along shocks, and supernova bubbles. ISM turbu- 
lence is a major driver of large-scale star formation: tur- 
bulent compression can incre ase the gas density and initi- 
ate star- forming instabilities (|Elmegreenlll993l . 120021 ). This 
can happen in the spiral arms of modern disk galax- 
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ies (jKlessenl |200l|), and there is recent evidence that 
very active star formation relates to increased turbulence 
and associated instabilities i n both primordial disk galax- 
ies (lElmegreen et al.l 120051: iForster Schreiber et al.l 120061 : 
iBournaud et all 120071: iDekel et all I2009T) and st arbursting 



mergers IjBournaud et a l. 2008a; 



Tevssierl l20ld h On the 



other hand, gas turbulence can also disrupt gaseous struc- 
tures faster than they gravitationally collapse, which could 
potentially stabilize gas dis ks and quench the star formation 
activ ity of some galaxies (jMartig et al.l |2009| ; IDekel et al.l 
l2009t ). 

Nevertheless, the large-scale properties of ISM turbu- 
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lence are still poorly understood and the ability of numeri- 
cal simulations to realistically model ISM substructures re- 
mains largely unexplored. Hydrodynamic simulations of ISM 
turbulence and f ragmentation in whole galaxy models have 
been performed (IWada et alJl200J; IWada fe Normanll200"7l; 



Taske r fc Brvanl 



20061 ; iTasker fc Tanl 120091 ; lAgertz et all 



2009a) but have not been compared to observational 



straints yet, probably because of lacking the required level 
of realism, both in spatial resolution and in global dynami- 
cal modeling (no old stellar components and live dark halo 
were considered), leading to a restricted dynamical range 
of structure sizes. All these models are characterized by a 
probability distribution function (PDF) of the gas density 
that has a log-normal shape, resulting from the combina- 
tion of various structures including density waves, dense gas 
clouds, supernovae bubbles. The levels of turbulent speed in 
these models seem realistic for both nearby and high-reds hift 
galaxies (see lAgertz et al.l l2009al ; ICeverino et al.l I2OIO1 . re- 
spectively). Whether a realistic spectrum of structures of 
various scales is reproduced in modern hydrodynamic mod- 
els remains however unknown. 

Observations suggest that turbulent motions are initi- 
ated by gravitational instabilities on a relatively large scale 
(the Jeans scale, which is 1/fc ~ 100 pc for wavenumber k 
in nearby disk galaxies, Elmegreen et al. 2003) and that a 
turbulent cascade develops towards smaller scales. This tur- 
bulence should be three-dimensional since the Jeans scale 
length is expected to set the gas disk scale height. Turbulent 
motion s could be tr iggered by gravity even inside molecular 
clouds (JFieldet al.ll2009l ). 

Observations of the power spectra of ISM emission from 
nearby galaxies are consistent with this. On scales smaller 
than ~ 100 pc, the power spectrum is steep, ~ —3, as it is 
for 3D Kolomogorov turbulence, while on scales larger than 
~ 100 pc, the power spectrum is flatter by about 1 in the 
slope, with a power-law exponent of ~ —2. This has been 
observed in the Large Magellanic Cloud (LMC) (Elmegreen, 
Kim, & Staveley-Smith 2001; Block et al. 2010), and in NGC 
1058 (Dutta et al. 2009). This tentatively suggests a two- 
dimentional regime for turbulent motions on large scales, 
because the break in the density distribution occurs at the 
expected disk scale height. The nature of the large-scale mo- 
tions is uncertain, however, because the three-dimensional 
velocity field in individual galaxies is not observable. In the 
Small Magellanic Cloud, the power spectrum is steep every- 
where (Stanimirovic et al. 2000), presumably because there 
is no thin disk on the line of sight. 

In this paper, we explore the properties of ISM turbu- 
lence in simulations of isolated disk galaxies, performed with 
an adaptive grid hydrodynamic code (AMR) with a maximal 
resolution of 0.8 pc resolution. We use a "pseudo-cooling" 
approach to model the equilibrium of gas between various 
heating and cooling processes, based on a barotropic equa- 
tion of state (EoS), and perform models with and without 
supernovae feedback. Our main model has a mass distribu- 
tion and rotation curve representative for the LMC disk, so 
as to compare with observations from Block et al. (2010). 
However we do not try to reproduce the detailed individ- 
ual structures of the real LMC, and we also analyze models 
with higher disk masses and gas fractions. We show that a 
large range of substructures spontaneously form in the ISM, 
with a power spectrum quite consistent with observations 



of the LMC from very large scales (> 1 kpc) down to very 
small scales (< 10 pc). A break in the power-law density 
power spectrum is observed at the gas disk scale height. 
Studying the velocity structure on various scales, we show 
that this corresponds to a transition from a 3D regime to a 
(quasi-) 2D regime in the gas motion, as initially proposed 
from observations (Elmegreen et al. 2001). Even if large- 
scale forcings are probably important, the two-dimensional 
regime on large scales cannot be interpreted just as density 
waves and streaming motions in a rotating disk, and evi- 
dence for an inverse cascade of turbulence in a quasi-2D disk 
is proposed. The same properties are recovered in a model 
of a massive high-redshift disk galaxy, which has a high gas 
fraction (50%) and strong turbulence (V/a ~ 4), suggesting 
these properties are relatively universal for gaseous galactic 
disks. 

Analyzing the vertical structure of the disk, compared 
to the properties of gas turbulence, we suggest that the tur- 
bulent motions at all scales are mainly powered by gravita- 
tional instabilities around the Jeans scale length, which sets 
the disk scale height. Feedback processes such as supernovae 
explosions do not significantly change the initial statistical 
properties of the ISM. They are nevertheless fundamental 
to maintain a realistic ISM structure over a large number 
of disk rotations. Supernovae bubbles do not appear to set 
the disk scale height, but they expand up to the disk scale 
height that is maintained by gravitationally-driven turbu- 
lence: this breaks apart the densest self-gravitating clouds 
to maintain a steady state cloud population. 



2 SIMULATIONS 

2.1 Simulation technique 

We perform closed-box model of isolated galaxies, with phys- 
ical parameters described in se ction 12.21 We use the AMR 
code RAMSES (|Tevssierl [20021 ). Stars and dark matter are 
described with collisionless particles, evolved with a particle- 
mesh (PM) technique. The gaseous component is described 
on the same adaptive grid; hydrodynamic equations being 
solved with a second-order Godunov scheme. 

Our simulations use a box size of 26 kpc. The coarse 
level of the AMR grid, l m in = 9, corresponds to a (2 9 ) 3 = 
512 3 Cartesian grid, i.e. the minimal spatial resolution is 
(■min — 51 pc where e = W~ kpc is the cell size. The maximal 
refinement level in dense regions is set by l m ax — 15, i.e. 
an effective resolution equivalent to 32768 3 and a maximal 
spatial resolution of e m a X — 0.8 pc. As we want to resolve 
turbulent motions in the whole disk, we use a refinement 
strategy that ensures that the mid-plane of the initial gas 
disk is entirely resolved with I > 11 (e < 13pc), and with 
I — 12 (e = 6 pc) at the initial average surface density. The 
coarsest resolution levels are used only in the outer stellar 
bulge and halo. The resolution increases towards a cell size of 
0.8 pc in the inner disk and in dense substructures forming 
in the outer disk. Effective resolution maps are shown on 
Figure |31 The refinement is as follows: a cell of level I < lmax 
is refined if its gas mass is larger than m res =5x 10 3 Mq 
and/or the number of particles (stars and dark matter) is 
larger than n res = 8. 

Only the gas density is computed in the smallest cells 



© 0000 RAS, MNRAS 000, 000-000 



Net tooling rale 



■:<i 




" 



. 



-a o 

log n„ <cm' ! )i 



Figure 1. The "pseudo-cooling" Equation of State used in our 
simulations. Gas cooling and heating rates (in erg s cm -3 ) as a 
function of temperature and density is shown for 1/3 solar metal- 
licity, as computed with the CLOUDY code l lFerland et al.lll99a . 
see text). The EoS used in our simulations is shown by the dashed 
line. A low-density regime corresponds to the Virial temperature 
of hot gas in a diffuse halo. Denser gas follows an equilibrium given 
by the minimal cooling rate: an isothermal regime accounts for 
warm HI in the disc at 10 4 K, and denser phases corresponding to 
cool/cold atomic and molecular gas. The EoS follows the thermal 
equilibrium, and variations of metallicity and UV background do 
not largely shift this equilibrium, so the same EoS would remain 
sensible for higher metallicities or UV fluxes. A thermal support 
at very high density is added to this EoS to prevent artificial 
fragmentation (see text). 
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Figure 2. Rotation curve of our LMC-sized model, measured at 
t = 254 Myr. 



of the AMR grid with I — l ma x = 15, the mass density 
of particles being restricted to lmax, P art = 14 in the PM 
scheme. This ensure a gravitational softening of at least 2 pc 
for particles. Treating the particle density at the l ma x level 
could result in a low number of particles per cell, and induce 
two-body relaxation, in regions that are refined because of 
a high gas density without necessarily having a high density 
of stars and/or dark matter. Only a few l m ax cells at each 
instant would actually contain less than 3 particles per cell, 
as measured in our simulation results, but we prefer avoiding 
any spurious effect. 

Resolving the gas cooling and heating equation is nu- 
merically costly and would reduce the achievable resolution. 
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In order to resolve small-scale structures and gas turbulence 
at the best possible resolution, we model the cooling and 
heating processes using an EoS that fixes the gas temper- 
ature T as a function of its local atomic density n. Here 
we use an EoS modeling the equilibrium of gas between the 
main heating (UV heating and stellar feedback) and cool- 
ing ( atomic and dust cooling) processes (see also iTevssieil 
120101 ) . This EoS is a fit to the average equilibrium tempera- 
ture of gas at one third solar metallicity: for densities 10~ 3 < 
n < 0.3 cm" 3 , T(n) = 10 4 K. For densities n > 0.3 cm" 3 , 
T(n) = 10 4 (n/0.3)" 1/2 K. Densities n < 10" 3 cm" 3 corre- 
spond to gas in gravo-thermal equilibrium in the halo, at 
the Virial temperature T(n) = 4 x 10 6 (n/10" 3 ) 2/3 K. 

Gravitational instabilities are a major driver of turbu- 
lence in galactic disks, both by directly triggering gas mo- 
tions and by fueling star formation that can also pump tur- 
bulent energy thr ough feedback pr ocesses. This is known 
from observations (Elmcgrccn 2002J), and models similar to 
ours but scaled to massive spiral galaxie s IjWada et al. 2002; 
iTasker fc Tar3l2009tlA~gertz et al.ll2009al ). A fundamental as- 
pect in our models is thus to avoid numerical instabilities 
and artificial fragmentation: a gas medium that should the- 
oretically be stable, tends to be unstable in numerical mod- 
els when the Jeans length is described by a small number of 
spatial resolution elements (jTruelove et al.lll997l ). It is gen- 
erally admitted that this effect is avoi ded if the Jeans length 
is resolved by at lea st a few elements (jMachacek et alj|2003i : 
lAgertz et a l. 2009a, e.g.). We then force the grid refinement 
to ensure that the thermal Jeans length of the gas is always 
resolved by at least 4 cells - and hence the actual ther- 
mal+turbulent Jeans length is resolved by an even larger 
number of cells). At the l ma x grid level, we impose a temper- 
ature threshold computed to keep the thermal Jeans length 
larger than 4 x e max . This temperature floor, added to the 
EoS at very high density (~ 10 5 Mq pc" 3 and above) can be 
considered as a sub-grid model for the unresolved turbulent 
motions at scales smaller than the maximal grid resolution. 
We can then consider that artificial fragmentation is gen- 
erally avoided: simulations of galactic disk fragmentation 
imposing higher numbers of resolution elements per Jeans 
length do not show significant differ ences in the properti es 
of gas stability and turbulence (e.g.. lCeverino et al.ll2010l ). 

One of our models includes only disk self-gravity, but 
the other model also includes star formation and feedback. 
A threshold for star formation is implemented at a gas mass 
density p = 5000 atoms cm" 3 ~ 115 Mq pc" 3 . This choice 
ensures that stars form only in the densest gas structures 
formed in our models but that the density PDF remains 
correctly sampled around and above this threshold (see Fig- 
ure [SJ. Above the threshold, the specific star formation rate 
(SFR per gas mass unit) is defined as the product of the 

efficiency and the local free-fall time* 



3tt 
32Gp' 



We use an ef- 



ficiency of 10%, which was chosen to give a gas consump- 
tion timescale of 2 Gyr on the results of the first simulation 
without star formation. Note that lower resolution models 
usually need to combine a lower efficiency with a lower den- 
sity threshold for star formation (see discussion in Teyssier 
et al. 2010). 

We use the kinetic feedback model described by 
JDubois fc Tevssierllirjoil ). A 20% fraction of the mass of 
stars formed at a given timestep is assumed to be in stars 
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that will explode into supernovae. 15% of the typical energy 
of each supernova, 10 51 erg, is re-injected into the gas, as- 
suming that the remaining energy is radiated away. The en- 
ergy is re-injected in the form of radial velocity kicks around 
the supernova, in a gas bubble of radius 3 pc (see Dubois & 
Teyssier 2008 for details). 

We use outflow boundary conditions for the fluid, and 
isolated boundary conditions for gravity 

2.2 Model Parameters 

Our simulations are not aimed at reproducing the exact mor- 
phology of the LMC, but at studying the properties of tur- 
bulence in a galactic disk with general properties similar to 
the LMC. To this aim, we initialize an exponential stellar 
disk containing 4x 10 6 particles, of total mass 3x 10 9 Mq ra- 
dial scale-length of 1.5 kpc, and truncation radius of 3 kpc. 
A non-rotating bulge of 3 x 10 8 Mq with radius 500 pc, 
made-up of 4 x 10° particles, is added. The stellar velocity 
dispersions initially correspond to a stellar Toomre parame- 
ter of Qs = 1.5, so that stars are stable against axisymmetric 
instabilities but unstable to bar forming instabilities - a bar 
and a pair of spiral arms form spontaneously in our model. 

The dark matter halo contains 4 x 10 6 particles, has 
a total mass of 5 x 10 9 Mq and follows a Burkert profile 
of core radius 3.5 kpc, truncated at 10 kpc. The initial gas 
disk is exponential, with a scale-length of 3 kpc and trun- 
cation radius of 5 kpc, an exponential scale-height of 50 pc 
truncated at 500 pc, and a gas mass of 6 x 10 8 Mq. The 
gas disk is initially purely rotating, with no macroscopic ve- 
locity dispersion. The initial densities in the disk lie within 
the isothermal part of our EoS, so a thermal support of 
~ 10 km s" 1 (T = 10 4 K) is present all over the initial disk, 
which can be considered as a model for the turbulent support 
of typically a few km s _1 in real disks. This support ensures 
an initial Toomre parameter for the gas Qo — 1 — 1.5. Once 
dense structures form in the course of the simulation, the 
temperature decreases (following the EoS) and the thermal 
support is gradually replaced by a turbulent support. Out- 
side the initial disk truncation, the AMR grid is initialized 
with a background density of 10~ 6 atom per cm 3 . 

The rotation curve in this model is shown on Figure [21 
It is broadly similar to the LMC rotation curve (Kim et 
al. 1998), with almost solid body rotation up to radii of at 
least 2 kpc, and gradually flattening with a circular velocity 
~ 70 km s _1 at large radii. We study the internal dynamics 
of the LMC; the interaction with the Milky Way and its 
hot gas halo is not modelled. Ram pressure was found to 
compress the leading edge of the HI and Ha disk, with little 
effect on the inner star-forming disk (Mastropietro et al. 
2009; see also Tonnesen & Bryan 2009 on the effect of ram 
pressure on a cloudy ISM). 




Figure 3. Top: Gas density maps of the LMC-sized gas disk 
with feedback, at t = 254 Myr, in a face-on projection (7x4 kpc 
snapshot) and a side-on projection (7 X 0.8 kpc snaphot). Bot- 
tom: Map of the AMR grid refinement level in the mid-plane, 
with I = 15 in white and I = 11 in dark blue. 



3 RESULTS 

The main simulation analyzed in Sections 3.1 and 3.2 is the 
LMC-sized model with stellar feedback. We later compare 
to the simulation without star formation and feedback in 
Section 3.3, to better highlight the role of feedback in the 
global structure of the ISM. 



3.1 ISM structure and density power spectrum 

We measured the velocity dispersion a in the gas disk, and 
found that both the mass-weighted average value < a > 
and the standard deviation Act rapidly increase in the first 
200 Myr and do not significantly evolve after t ~ 200 Myr, 
indicating that a steady state in the turbulent regime has 
been reached. We thus mostly analyze our models around 
i=250 Myr. 

Spectral analysis of the gas distribution was performed 
on a sub-box of size 13 kpc (the outer regions of the full 
simulation box contain only the outer regions of the dark 
matter halo). Throughout the paper, the wavenumber k = 1 
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Figure 4. Zoomed-in face-on and side-on density maps of the LMC-sized model with stellar feedback at t = 254 Myr. The horizontal size 
of these figures is 3 kpc. The darkest regions visible on these projections have surface densities of 35 Mq pc~ 2 : a low density fountain 
above the disk plane is not shown on the edge-on projected. 



is for a wavelength of 13 kpc, which results in a wavenumber 
unit u k - 4.70 x 1CT 4 pc" 1 . 

Figure [3] shows the gas density distribution for this sim- 
ulation at t — 254 Myr. The response to a stellar bar, and 
some long spiral arms are seen, together with more chaotic 
structures such as shorter arms, shocks, dense clouds, bub- 
bles, etc (Fig. [4]). We fit a sech 2 vertical profile to the gas at 
radii smaller than 5 kpc, integrated over the plane, and find 
a scale height of 207 pc. The scale height is relatively con- 
stant, with moderate flaring, and varies from 187 pc inside 
the central kpc to 226 pc for 3 < r < 5 kpc. 

The density PDF, integrated over the central disk with 
a radius of 5 kpc and height of 2 kpc, is shown on Figure [5] 
This PDF has a log-normal shape, truncated at quite high 
density (~ 10 6 cm -3 ) which is enabled by our high spatial 
and mass resolution. If we write the probability distribution 



for the density 



P{p) 



-o. 



•>*= 



oc e a 5 , 

then the Gaussian dispersion is A ~ 2.92. 

The power spectrum of the face-on gas column density 
is shown on Figure [6] at the same instant and later ones. 
No significant time evolution is found. The power spectrum 
shows a double power-law with a break at wavelengths of 
about 150 pc, which is about the measured gas disk scale- 
height. The power law for large structures has a slope of 
~ —1.9 and the slope for small structures is ~ —3.1. The 
slope steepens for the smallest scales below 5-10 pc, which 
is discussed later in section 3.3. 

The simulation without star formation and feedback 
shows a globally similar density power spectrum (see 
Sect. 3.3 and Fig. 14) with a double power-law shape with 
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Figure 5. Probability distribution function (PDF) of the gas 
density p at t = 254 Myr in the model with feedback. The dashed 
curve is a log normal profile centered on p ma x = 53 cm -3 and 
of standard deviation A=1.27 in a base-10 log PDF (see text for 
details). The dotted line is for the simulation with reduced reso- 
lution (section 13. 5 f . This PDF was measured within a cylindrical 
box of radius 5 kpc and height 2 kpc. The peak at low densities 
corresponds to low density gas in the hot halo. 
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Figure 6. Power spectrum of the face-on gas surface density 
in the model with feedback, at t = 254 Myr (dark), 261 Myr 
(green), 268 Myr (blue) and 275 Myr (red). The wavenumber 
unit is 4.70 x 10 -4 pc -1 . 



slopes of — 2.9 on small scales and —1.8 on large scales. The 
transition occurring around the gas disk scale-height. Only 
the smallest scales below 10 pc significantly differ from the 
model with star formation and feedback (see section 3.3). 

Our simulations produce a density structure that has 
about the power spectrum observed for the ISM in the 
LMC and other nearby galaxies. The power spectrum shows 
a break between a —2 power law and a —3 power law, 
with the transition occurring at scales around the gas disk 
scale- height, as usually speculated in observations (e.g., 
lElmegreen et al.ll200ll ; iDutta et al.ll2009r i. In the following, 
we study the velocity structure to understand the nature of 
the motions in the two regimes of the density distribution. 
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Figure 7. Velocity power spectrum for the three separate com- 
ponents V r , Vq and V z , at t = 254 Myr (LMC-sized model with 
stellar feedback) . The conversion of wavenumbers into linear size 
is as on Fig. [6] 




Figure 9. Face-on gas density snapshot at t = 268 Myr, with 
a 7 X 4 kpc field of view, in the LMC-sized model with stellar 
feedback. 



3.2 Velocity structure 

3.2.1 Velocity fields and power spectra 

Maps of the in-plane radial velocity component V r and per- 
pendicular component V z are shown on Figure[8]for the same 
snapshot as Figures [3] and [4] all at t — 254 Myr. Another 
instant (i = 268 Myr) is shown on Figure [5] for the gas den- 
sity and Figure [10] for the velocity maps. The velocity maps 
are shown for the total velocity components, and for the low 
wavenumbers k < 50 (A > 300 pc) and the high wavenum- 
bers k > 200 (A < 70 pc) separately. This separates the two 
regimes identified on the density power spectrum, where the 
transition occurs at k ~ 100 or A ~ 150 pc. These maps were 
built by zeroing the low-fe or high-fe components in the result 
of a Fourier transform of the velocity field, and recovering 
the velocities through an inverse Fourier transform. 

The power spectrum of the three velocity components 
V r , Ve and V z is shown on Figure [7] A single power law is 
found for the in-plane components. The perpendicular ve- 
locity V z follows the same power law on small scales, but its 
power spectrum flattens on large scales, with a transition 
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Figure 8. Gas velocity fields at t = 254 Myr. The in-plane V r component is shown in the left column, and the perpendicular component 
V z us shown to the right. Both are mass-weighted averaged along the line of sight. We show the total velocities (top), and the large-scale 
(k < 50, middle) and small-scale (fc > 200, bottom) components, after zeroing the other wavenumbers in a Fourier decomposition. The 
break in the density power spectrum is at k ~ 100 and separates the two components shown separately here. The same colorbar is used 
for all maps. The field correspond to the face-on gas density image shown on Figure [3] in 7x4 kpc snapshots (LMC-sized model with 
stellar feedback). 



at about the same scale length 
structure spectrum. 



^150 pc) as in the density 



3.2.2 A 3D cascade of turbulence on small scales 



These differences in velocity power spectra correspond 
to higher power for the long-wavelength in-plane compo- 
nents than the long-wavelength perpendicular components 
of velocity. There are large in-plane motions at long wave- 
lengths from spiral and bar-driven gas flows and relatively 
little perpendicular motions at long wavelengths to accom- 
pany them. 



On scales smaller than the gas disk scale height (100-200 pc), 
the ISM density has the density power spectrum expected 
for a 3D cascade of turbulence. The three-dimensional na- 
ture of the associated motions is confirmed by the high- 
wavenumber velocity maps, which show about the same 
amplitude for the in-plane motions and the perpendicular 
one. Typical average values of < V z /V r > on various scales 
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v z 




total 



k < 50 



k > 200 




Velocities (km s *) 
Figure 10. Same as Figure[8]but at t = 268 Myr, corresponding to the Figure l9l density snapshot. 



are given in Table Q] The small-scale motions are globally 
isotropiqjo 

These results are overall consistent with a 3D cascade 
regime on the small-scale part of the density spectrum, 
initiated at scales around 150 pc, which is the gas disk 
scale height and the average Jeans length. This is natu- 
rally the case if the most gravitationally unstable scale, 
namely the Jeans length, corresponds to the disk thick- 



1 The slightly higher values for the vertical velocity component 
correspond to low-density gas flows ejected above the disk plane 
by supernovae winds; the simulation without feedback has values 
of < V z > even closer to < V r > for high wavenumbers. 

2 Note that V r is slightly larger than Vg , so the in-plane motions 
are not perfectly isotropic. This is expected for instance if the gas 
tends to follow epicyclic motions, for which the ratio of the radial 
to tangential velocity dispersions is 2Q/k. 



ness. Indeed, the gas disk scale height is expected to be 
set by the gaseous Jeans length in near by disk galaxies 
iJElmeereen et al.ll20oJ ; iDutta et aP 12009? ) : this is also the 
case in high-redshift disk galaxies, which are more turbu- 
lent, have large r Jeans lengths, and have thick e r gas disks at 
the same time C Elmegreen & Elmegrcen 200 a: iGenzel et al.l 



l200d ; iForster Schreiber et alJl200q ; lBournaud et al.ll2008bh" 
We checked this in our simulations by computing a mass- 
weighted average of the velocity dispersion a av (8 km s _1 ), 
the average surface density S al) , and then an "effective" 

Jeans length l/fcj ea ns,cft = ^ G ^ = 213 pc (locally, the 
Jeans length can significantly differ from this average value, 
as the density and velocity dispersion vary). This is consis- 
tent with the disk scale height being set by the Jeans length. 
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Table 1. Mass-weighted average values of < V z /V r > for the 
large scale, low scale, and total velocities (LMC-sized model with 
feedback) . 

Time 

small scales (k > 200) 
large scales (k < 50) 
total velocities (all k) 



254 Myr 


268 Myr 


1.13 


1.07 


0.16 


0.14 


0.92 


0.89 




Figure 11. Vorticity map for the large-scale (k < 50) com- 
ponents, at t = 254 Myr (LMC-sized model with stellar feed- 
back). The scale ranges from — 10 10 pc 2 Myr -1 (black) to 
10 10 pc 2 Myr -1 (white) with the same color bar as velocity fields. 



3.2.3 A 2D inverse cascade on large scales? 

The low-wavenumber (large-scale) structures have a flatter 
density power spectrum. They also have different kinematic 
properties. The large-scale component of the velocity fields is 
dominated by in-plane motions. Perpendicular motions are 
much weaker. The low-wavenumber component of V z shows 
only modest peaks in Figure[5] and these peaks are generally 
found in low-density regions (they generally correspond to 
gaseous fountains above the disk plan, rather than vertical 
motions in the mid- plane). 
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Figure 12. Enstrophy spectrum at t = 254 Myr for the LMC- 
sized model with stellar feedback. The wavenumbers are divided 
by a factor of 2 compared to the density and velocity power spec- 
trum (Fig. [6]and[7]l. Scales larger than 100-200 pc (here k ~ 2) 
have a -1 power law spectrum: the amplitude of vortices of various 
sizes is just as expected for two-dimensional turbulence. 



The motions corresponding to the large-scale regime of 
the density power spectrum are quasi-2D, highly anisotropic 
(Table [T]). This is not because the disk rotation dominates 
these motions: this applies to V r as much as Vg, and the 
< V z /V r > ratio is almost unchanged when we remove 
the lowest wavenumbers k ^ 3 tracing global rotation. The 
power spectra of the three velocity components (Fig. [7]) ac- 
tually shows that the perpendicular velocity is much lower 
than the in-plane components for all wavelengths larger than 
the disk scale-height. 

The ~ —2 slope of the large-scale branch of the density 
power spectrum and the quasi-2D nature of the associated 
motions suggest that we are observing something like a 2D 
inverse cascade of turbulence (Kraichnan 1967). We never- 
theless have to explore the nature of these motions more 
accurately, to ensure that these properties are not a con- 
spiracy from density waves and associated gas flows, which 
would then not be turbulent motions. 

The large-scale components of the in-plane motions 
(Figs [8] and I10[) show a central m — 2 mode, which is 
the characterist ic response to a bar+arms system of den - 
sity waves (e.g.. lCombes fc Gerirj|l985l : lAthanassoulall2000h . 
This response is however quite asymmetric, and does not 
dominate the large-scale motions outside of the central kpc. 
Overall, outside the central kpc, chaotic motions dominate 
the large-scale velocity components (of course superimposed 
on the disk rotation pattern). In particular the comparison 
of the gas density and velocity maps shows large vortices 
around dense gas clumps: for instance, the densest clump 
seen to the left of Figure [3] is clearly associated with an 
in-plane eddy, with positive and negative peaks of V r sur- 
rounding the position of this dense clump. Such signatures 
are reminiscent of a Rosby W ave instability (|Lovelace et al.l 
ll999l : IVarniere fc Taggerll20061 ). which is obser ved in pure 2D 
disks but also in thick disks or 3D systems IjMeheut et al.l 
120101 ) and could thus take place in our quasi-2D system. 
Such motions would be quasi- 2D turbulence, which is more 
than a large-scale flow in a density wave, but which may be 
induced by a density wave of other large-scale forcings. 

The vorticity map of the large-scale velocity component 
(k < 50), after subtraction of the m = rotation pattern, 
is shown on Figure [TT] It shows a number of vortices with 
various sizes and strengths throughout the disk. The power 
spectrum of the enstrophy is shown in Figure [T2j Enstrophy 
is the integral of the square of the vorticity: £ — J 1 1 V x v 1 1 2 . 
The enstrophy has a power-law power spectrum through- 
out the entire range of 2D scales, like the in-plane motions 
shown before. The slope of this power spectrum is -1, just as 
expected for 2D turbulence, as found in two -dimensional ex- 
periments (e.g.. |Petersen et al.ll2006l . 120071 ). This enstrophy 
spectrum means that the quasi-2D part of the density power 
spectrum is made-up of vortices with size and strength dis- 
tributed as expected for 2D-turbulence: this is unlikely to 
simply result from non-turbulent flows along density waves. 

The properties of the gaseous motions on scales larger 
than the disk scale height thus contain a component from 
turbulence in a quasi-2D medium. The origin of this tur- 
bulence could be a combination of long-range forcing from 
spiral waves and other global disturbances, and an inverse 
cascade from smaller scales, where energy is put into the 
ISM by gravitational instabilities on the Jeans length and 
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stellar feedback. Energy injection at the Jeans length is con- 
sistent with our observation of vorticity around the gravi- 
tating clumps and with the quasi-2 D power sp e ctrum on 
scales larger than the Jeans length. IWada et al.l (|2002T ) al- 
ready mentioned that the properties of turbulent motions in 
their galaxy models were consistent with an inverse cascade 
of turbulence. Their models were purely two-dimensional, 
while we here suggest a similar inverse cascade in a three- 
dimensional gas disk, which is thin but has a well resolved 
scale height. 

The large-scale motions probably cannot simply result 
from a single inverse cascade with energy input only at the 
Jeans scale. Indeed, the growth rate of gravitational pertur- 
bations in a disk decreases only as 1/vA where A is the scale 
of the perturbations, so that long-range forcing at scales 
larger than the typical Jeans scale length remains significant. 
At least, the gravitational coupling of stel lar and gaseous 
comp onents and the bi-component stability l|Jog fc Solomon! 
Il984f ) should result in a large range of unstable scales, rang- 
ing ty pically from the g aseous Jeans scale to the stellar Jeans 
scale (lElmegreenlll995T I. There should then be a range of un- 



stable scales injecting energy into the inverse turbulence cas- 
cade. While the expected properties of 2D turbulence, such 
as the enstrophy spectrum, are recovered in our simulations, 
one can note that the slope of the density power spectrum 
in this regime is slightly higher than expected from pure 
two-dimensional turbulence (a —2 power law, instead of the 
—5/3 exponent expected for purely two-dimensional and un- 
compressible turbulence). The non-zero thickness of the gas 
disk and the relative importance of long-range forcing in a 
gaseous+stellar disk might cause this. 

The role of long-range gravitational forcing and the 
coupling of the gaseous and stellar components also ap- 
pear by comparing our sim ulations with the AMR mod- 
els in fTasker fe Tanl (|2009h . The stellar disk is modeled 
with a rigid analytical potential in their simulations. Gas 
clumps form everywhere in their disk with a relatively 
peaked size distribution and a constant separation of about 
one Jeans scale. The Jeans length is also a characteristic 
scale in our model: it sets the disk scale height and gas 
clouds along spiral arms are also often separated by Lj 
(Fig. Q, but gaseous structures overall form over a wider 
range of sizes in our models. The inclusion of both stellar 
and gaseous gravity on large scales appears crucial to re- 
alistically reproduce the full spatial range of the main ISM 
structures from long spiral arms and large vortices to smaller 
clumps and shocks, even if ISM structures on small scales 
can be realistically reproduced w ithout stellar gravity (e.g., 
Ide Avillez fc Breitschwerdtll200^ . 



We have analyzed above the LMC-sized simulation with 
star formation and stellar feedback. In this simulation, the 
properties of ISM turbulence seem mainly driven by gaseous 
and stellar gravity. However, the comparison to the model 
without stellar feedback below will highlight a fundamental 
role of feedback in maintaining a steady state in the density 
distribution of the ISM. 



3.3 Role of star formation feedback in gas disk 
properties 

The "gravity-only" simulation (without star formation and 
feedback) shows a relatively similar distribution of gaseous 
structures, and a relatively similar density power spectrum 
compared to the model with feedback (see Figures [13] and 
I14[) , at least at t — 238 Myr and at scales larger than 5-10 pc. 
Besides suggesting that our earlier results do not crucially 
depend on a specific feedback scheme, it also indicates that 
the ISM structuring and the pumping of turbulent energy 
into the turbulent cascades results mainly from gravitational 
processes and is not primarily driven by supernovae explo- 
sions. The disk scale-height and the disk substructures are 
unchanged in this gravity-only model. 

Feedback processes nevertheless play a fundamental role 
in maintaining this gravity-driven turbulence distribution. 
We can indeed note already at t — 238 Myr in Figure [Hl that 
the ISM power spectrum shows an excess of very small-scale 
structures (below 5-10 pc) and the gas density map shows 
the associated small and dense gas clumps. Over longer 
timescales the power spectrum becomes less realistic as the 
large scales are gradually emptied and the bump at very 
small scale increases: gas dissipates its energy and accumu- 
lates in tiny dense bullets from where it is never removed 

(Fig. El). 

The role of feedback appears when comparing the den- 
sity power spectra in our two models: feedback disrupts 
the dense smallest-scale structures, and supernovae bub- 
bles expand up to the disk scale height. Density maps show 
supernovae-driven bubbles, all of which are smaller than 
or comparable to the disk scale-height. But this does not 
mean that the size of supernovae bubbles regulates the scale 
height: the measured scale-height and the break in the den- 
sity power spectrum are not changed by feedback; the scale 
height is mostly regulated by gravitational processes, and 
is about the Jeans length. But feedback disrupts structures 
on the smallest scales and expands them up to the Jeans 
length, where gravitational processes can initiate a new cy- 
cle of 3D turbulent cascade towards small scales and 2D 
inverse cascade towards large scales. 

Overall, the gas disk structure and its scale height are 
largely regulated by gravitational instabilities, but gas struc- 
tures and energy dissipation on the smallest scales are also 
regulated by star formation feedback. The break in the 
face-on density power spectrum and the scale height mea- 
sured on edge-on projections both correspond to the average 
Jeans length. 2D and 3D motions are initiated by instabili- 
ties around this scale length and large-scale forcings. Feed- 
back processes prevent the otherwise inevitable accumula- 
tion of gas in tiny and dense bullets, disrupt small-scale 
structures and refill the turbulent cascades initiated at the 
Jeans length. This cycle maintains the ISM density distri- 
bution in a steady state. 

Observations have already suggested that turbulent 
motions are primarily driven by gravitati onal instabilities 
such as spiral arms and molecular clumps (|Elmegreen et al.l 
120031 ). Based on the analysis of the density distribution, 
these authors pointed out the role of star formation feed- 
back in brea king apart sm a ll dens e structures to preserve a 
steady state. lAgertz et al.l (I2009al ) also proposed from sim- 
ulations that the main driver of ISM turbulence could be 
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Figure 13. "Gravity-only" LMC-sized model without star forma- 
tion and feedback: face-on view of the gas density at t = 343 Myr. 
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Figure 15. Gas surface density power spectrum in a massive 
gas-rich galaxy model. 
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Figure 14. Face-on gas density power spectrum of the "gravity- 
only" model without stellar feedback at t = 258 (solid lines) and 
343 Myr (dashed, shown at small wavelengths only for clarity). 



gravitational forcing, with nevertheless a significant contri- 
bution from star formation (see also Dib, Bell & Burkert 
2006). 

Gravity-only models without star formation and feed- 
back initially produce a realistic density distribution in the 
ISM over a couple of disk rotation periods. But on longer 
timescales and over large numbers of disk rotations, the in- 
clusion of feedback in numerical models appears necessary to 
maintain a realistic density distribution. However, this ap- 
plies only to high-resolution simulations where gas cooling 
is modeled down to ~ 100 K. Lower-resolution SPH simu- 
lations, and more generally models where the ISM is only 
modeled as a smooth stable gas at ~10 4 K or more (e.g., the 
stabilizing EoS by Springel & Hernquist 2003 and Robert- 
son et al. 2004), would not be affected by catastrophic gas 
dissipation if they do not include feedback sources, as they 
do not resolve the turbulent and unstable nature of the star- 
forming ISM. 



3.4 Extension to massive and gas-rich galaxies 

The results described so far were for a low-mass disk galaxy, 
scaled to the LMC properties. To explore whether the high- 
lighted properties can extend to more massive disk galaxies 
and/or disk galaxies with different gas fractions, we have 
applied the same analysis to a high-red shift galaxy model 
described in lDelave et al.l (J2010 in prepl ). This is a massive 
galaxy (baryonic half-mass radius of 10 kpc and circular 
velocity of 280 km s _1 ) with a 50% gas fraction, modeled 
with a 2 pc resolution. As typical star-forming galaxies at 
z ~ 2, it has a high gas fraction, a very unstable disk with 
high turbulent speeds (a few tens of km s _1 ) and giant 
complexes of star formati on in a thick (h z ~ 1 kpc) ga s 
disk (see observ a tions by Elmegreen fc Elmegreenl {2006); 
lElmeereen et alj J2007I); iGenzel et all d2006T) and simula - 
tions bvlSemelin fc Combes! (l2005l);lBqurnaud et all (|2007h : 
lAgertz et all l|2009bl ); ICeverino et all (|201fj| )). ~ 

The density power spectrum is shown on Figure [151 It is 
consistent with a double power law, with the break occurring 
at wavelength around 1 kpc which is the disk scale height 
and the typical Jeans length in such galaxies. We measured a 
total ratio of velocity dispersions of < V z /V r >= 0.88, with 
< Vx/Vt >= 0.26 on large scales and < V z /V r >= 1.07 
on small scales. The large-scale motions are still signifi- 
cantly non-isotropic, even if the low wavenumber velocities 
are somewhat more isotropic than in our LMC-sized model, 
which could be caused by the disk being thicker. The prop- 
erties in the density distribution and chaotic motions found 
in our LMC-sized models thus seem to apply, at least qual- 
itatively, to galaxies with different masses, rotation curves 
and gas fractions. 



3.5 Resolution requirements 

Our fiducial LMC-sized model with star formation and feed- 
back has been re-done with l m ax = 11, i.e. a spatial resolu- 
tion of 13 pc, which means that the disk scale height is re- 
solved with ~10 resolution elements, instead of 100-200 res- 
olution elements at the maximal refinement level l max = 15. 
The density threshold for star formation was reduced by a 
factor of 80 to keep the resulting star formation rate similar. 
The global density distribution of this reduced- 
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resolution model is not very different from the initial case, 
except that the highest density structures in the PDF are 
removed; otherwise the PDF is only moderately shifted to- 
wards lower densities and has about the same dispersion 
(Fig. [5J). However, the velocity structure is significantly dif- 
ferent. The initial model had almost isotropic velocity dis- 
persion a z /a r ~ 0.9 (the smallest scales begin the most 
isotropic ones). Interestingly, the reduced-resolution model 
has much less isotropic motions, with <r z /oy = 0.32 (at 
t = 262 Myr and over all wavenumbers; the ratio was ~0.9 in 
the high-resolution run). We thus note that a very high res- 
olution is required to resolve the isotropy or non-isotropy of 
turbulent motions in disk simulations, and that a few (~10) 
resolution elements per disk scale height are not sufficient. 
Our interpretation is as follows: vertical motions can be initi- 
ated or amplified when instabilities occur off the midplane at 
heights z>0 or z<0, when the Jeans length becomes locally 
smaller than the disk scale height (which happens in dense 
regions or low- velocity dispersion regions). This requires to 
resolve, locally, several Jeans length per scale height. Then, 
each Jeans length itself needs to be resolved by a least a 
few resolution elements: 4 resolution elements per thermal 
Jeans length is a strict minimum imposed in our model and 
our EoS generally results in 5-10 elements per thermal Jeans 
length, and the number of resolution per thermal+turbulent 
Jeans length would be even higher. Thus, capturing insta- 
bilities at z > and the associated triggering of vertical mo- 
tions can require a few tens of resolution elements across the 
gas disk scale height. When the disk thickness is marginally 
resolved, the non-isotropy of gas motions can be artificially 
increased at small scales. In this case, each gas clump is 
mostly spinning around its minor axis (see Agertz et al. 
2009a) generally aligned with the spin axis of the entire 
galaxy (but see Tasker & Tan 2009). 



4 SUMMARY AND CONCLUSIONS 

4.1 Properties of ISM turbulence, role of gravity 
and feedback 

This paper has presented an LMC-sized models for compar- 
ison with observations of the LMC (Block et al. 2010) and 
a high-redshift clumpy disk model. The LMC-sized model 
does not match all individual structures observed in the 
LMC - this is not the initial purpose and the LMC has 
many low-density ionised filaments, and other substructures 
in low-density regions, that are out of reach of our mass 
resolution. 

The model is successful in reproducing the double power 
law shape of the ISM density power spectrum observed in 
the LMC (Elmegreen et al. 2001 and Block et al. 2010) in a 
system with a similar mass, gas fraction and rotation curve. 
We analyzed the 3D velocity structure of the models, which 
could not be done in observations, to understand the origin 
of this density substructures. 

Our results indicate that a dual turbulence cascade 
could take place in our simulations, and by extension in the 
real ISM. The scale height of the gas disk is set by the Jeans 
scale length, which is the most gravitationally unstable scale. 
Gravitational instabilities at/around this scale length inject 
energy in turbulent motions. These motions follow a classical 



3D cascade of isotropic turbulence towards smaller scales, 
on which they form a hierarchy of gas clouds, GMCs and 
substructures therein. Turbulent motions on scales down to 
and inside GMCs thus appear to be from self-gravity. Nev- 
ertheless, stellar feedback plays a major role in disrupting 
dense structures on the smallest scales and in maintaining 
the turbulence cascade in a steady state. 

On scales larger than the gas disk scale height, very 
anisotropic turbulent motions with a shallower density 
power spectrum are observed. Long-range gravitational 
forces are probably not negligible, but the vorticity and en- 
strophy properties on these large scales are as expected for 
(quasi-) 2D turbulence. The large scale motions are not lim- 
ited to global disk rotation and streaming flows along density 
waves. This inverse cascade appears to be initiated mainly 
by the same gravitational instabilities at the Jeans length as 
the direct 3D cascade, because it is present with and without 
stellar feedback. 

While models without self-gravity do not find a char- 
acteristic scale for energ y injection into ISM turbulence 
l|Joung fc Mac Low 2006), we identify a major role of grav- 
itational processes around the Jeans length. Large-scale 
forcing and stellar feedback at small scales are important 
processes in regulating the ISM properties too, but the 
Jeans length is the scale at which chaotic gas motions 
change from quasi-2D to 3D isotropic turbulence, and at 
which a break occurs in the density power spectrum. Log- 
normal gas density PDFs have been found in models of 
supernovae-driven ISM turbu lence without self-gravity (e.g. 
Ide Avillez fc Mac Low! 120021 ). but without having the ob- 
served power spectrum reproduced. The power spectrum 
analysis in our models and in LMC observations rather sug- 
gests that the energy injection in the ISM is largely from 
gravitational processes (such as gravitational instabilities 
and inward mass accretion , lElmegreen fc Burkertl [2010J; 
iKlessen fc HennebehafeOlOT ) with a regulating role of feed- 
back on small scales. 

We have shown that similar properties for ISM sub- 
structures and turbulence, at least qualitatively, are found 
in models of galaxies with higher masses and higher gas 
fractions. This includes the extreme case of high-redshift 
(z ~ 2) galaxies with thick, highly turbulent and very 
clumpy gaseous disks. An interesting implication for ob- 
servers is that the isotropy of the gas velocity dispersion is 
strongly dependent on the scale, and hence on the resolution 
of an observation. 



4.2 Modeling ISM substructures in galaxy 
simulations 

We have explored ISM modeling with a barotropic equation 
of state (EoS) corresponding to gas in thermal equilibrium 
between a standard radiation field with the main atomic 
and fine-structure cooling processes. We have shown that 
this EoS employed in an AMR hydrodynamic code with a 
quasi-Lagrangian refinement strategy reproduces a realistic 
power spectrum for the ISM density distribution. It does not 
form a truly multiphase medium, since the pressure is only 
a function of density. It nevertheless a clumpy ISM, with 
cold and dense clouds embedded in a warmer diffuse phase, 
and the resulting density distribution is consistent with that 
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of the real, multiphase ISM. This technique is efficient to 
reach high spatial and mass resolutions: the computation of 
cooling and heating rates is not numerically costly, but the 
absence of out-of-equilibrium, colder or warmer gas in this 
model prevents the timescale to become extremely short, 
which would typically happen in parsec-scale models with 
complete cooling calculations. 

Hydrodynamic models with c ooling calculations pro- 
duce a log-normal PDF for the ISM dWada fc Normanll2007l; 
iTasker fc Tan]|2009l ; lAgertz et al.ll2009aT ). which is retrieved 
in our model over a larger dynamical range. The PDF 
alone does not guarantee that the correct mass fraction 
is distributed in structures of various sizes and densities. 
Our model with a mass and rotation curve similar to the 
LMC matches the density spectrum observed by Block et 
al. (2010) in the LMC. While all thermal processes are not 
explicitely computed, this at least means that the ISM mod- 
eled this way has realistic density distribution, with a correct 
number of gas structures of various sizes and a correct mass 
fraction therein, from the global galaxy-sized scales down to 
5-10 pc, i.e. molecular clouds and their main substructures. 

The ISM substructures are mostly powered by gravity- 
driven turbulence and instabilities around the Jeans scale 
length. Star formation feedback is however an essential in- 
gredient, not only by structuring the ISM into shells locally, 
but also by providing an extra source of energy is needed to 
disrupt the smallest and densest structures that tend to be 
steadily filled by gas dissipation into compact clumps. Grav- 
ity, hydrodynamics, and a realistic thermal model are ini- 
tially enough to produce a realistic structure distribution in 
the ISM: feedback seems relatively negligible over a couple of 
disk rotations, but becomes important over longer timescales 
to balance systematic gas dissipation. The source of energy 
that disrupt the structures at very small scales and refills the 
3D cascade and the inverse 2D cascade at the Jeans scale is, 
in our model is supernova feedback. Other sources such as 
HII regions, st ellar winds, and rad i ation pressure from young 
massive stars JMurrav et al.ll2010l : JKrumholz fc DekellboiOT ) 
could also contribute to maintain the ISM density distribu- 
tion in a steady state. The global density structure of the 
ISM remains relatively unchanged in simulations where the 
disk scale height is marginally resolved, but modelling accu- 
rately the quasi-isotropic turbulent motions on small scales 
is found to require a very large number (at least a few tens) 
of resolution elements per disk scale height. 
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